!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine block_interpolate (child, nx_c, ny_c, parent, nx_p, ny_p, start_x, start_y, ratio)
!==================================================================================================================================|

   implicit none

   integer(4), intent(in) :: nx_p, ny_p, nx_c, ny_c, start_x, start_y, ratio

   real(8), dimension(nx_p, ny_p), intent(in) :: parent
   real(8), dimension(nx_c, ny_c), intent(inout) :: child

   integer(4) :: i, j, i_p, j_p, ip1, im1, jp1, jm1
   real(8) :: lmt, eps, denominator1, denominator2

!----------------------------------------------------------------------------------------------------------------------------------|
   eps = 1.0d-12

   do j = 1, ny_c

      j_p = (start_y + j) / ratio
      jp1 = j_p + 1
      jm1 = j_p - 1

      do i = 1, nx_c

         i_p = (start_x + i) / ratio
         ip1 = i_p + 1
         im1 = i_p - 1

         denominator1 = parent(ip1, j_p) - parent(im1, j_p)
         if (denominator1 .eq. 0.0d0) denominator1 = eps

         denominator2 = parent(i_p, jp1) - parent(i_p, jm1)
         if (denominator2 .eq. 0.0d0) denominator2 = eps

         lmt = max ((parent(ip1, j_p) - parent(i_p, j_p)) * (parent(i_p, j_p) - parent(im1, j_p)), 0.0d0) /                        &
            denominator1 * 2.0d0 / ratio * (mod (start_x + i, ratio) - (ratio - 1) / 2.0d0) +                                      &
               max ((parent(i_p, jp1) - parent(i_p, j_p)) * (parent(i_p, j_p) - parent(i_p, jm1)), 0.0d0) /                        &
            denominator2 * 2.0d0 / ratio * (mod (start_y + j, ratio) - (ratio - 1) / 2.0d0)
         child(i, j) = parent(i_p, j_p) + lmt
      enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine block_interpolate
